

#ifdef OPENCL
ocl_t *setup_ocl() {

	ocl_t *ocl;
	ocl = calloc( 1, sizeof( ocl_t));

//	fprintf( stderr, "setup_ocl()\n");

	const char* const_kernel_code =
	"#define BLOCKSIZE	384	/* optimal on Fermi */                                                                            \n"
	"#define DAMP		2.1304f	/* exponential damping width */                                                               \n"
	"#define MATRIXSIZE	16                                                                                                    \n"
	"#define MAXFVALUE	1.0e14f                                                                                               \n"
	"                                                                                                                         \n"
	"__kernel void thole_iterative_ocl                                                                                        \n"
	"(                                                                                                                        \n"
	"  __global float4 *arg_ptr,                                                                                              \n"
	"  __global float  *basis,                                                                                                \n"
	"  __global float  *recip_basis,         																				  \n"
	"  __const  uint    N                                                                                                     \n"
	")                                                                                                                        \n"
	"{                                                                                                                        \n"
	"   __local float4 sposj[ BLOCKSIZE];                                                                                     \n"
	"   __local float4 smu[ BLOCKSIZE];                                                                                       \n"
	"	int bid, tid, bsize, gsize;                                                                                           \n"
	"        int i, j, k, mol_id;                                                                                             \n"
	"	__global float4 *pos, *estat, *mu_in, *mu_out, *eind_out;                                                             \n"
	"        float4 dr, dri, img;                                                                                             \n"
	"	float4 posi, posj, estati, mu, eind, mu_out_r = (float4)(0.0f,0.0f,0.0f,0.0f);                                        \n"
	"	float r, r2, r3, r5, ri, ri2;                                                                                         \n"
	"        float damp = DAMP, damp2, damp3, damping_term1, damping_term2;	/* exponential damping func */                    \n"
	"	float expr;                                                                                                           \n"
	"	float sw;							/* switch to avoid conditional */                                                 \n"
	"   float4 basis_r_a = (float4)( basis[ 0], basis[ 1], basis[ 2], 0.0f);                                                  \n"
	"   float4 basis_r_b = (float4)( basis[ 3], basis[ 4], basis[ 5], 0.0f);                                                  \n"
	"   float4 basis_r_c = (float4)( basis[ 6], basis[ 7], basis[ 8], 0.0f);                                                  \n"
	"   float4 recip_basis_r_a = (float4)( recip_basis[ 0], recip_basis[ 1], recip_basis[ 2], 0.0f);                          \n"
	"   float4 recip_basis_r_b = (float4)( recip_basis[ 3], recip_basis[ 4], recip_basis[ 5], 0.0f);                          \n"
	"   float4 recip_basis_r_c = (float4)( recip_basis[ 6], recip_basis[ 7], recip_basis[ 8], 0.0f);                          \n"
	"                                                                                                                         \n"
	"   /* dipole field tensor    */                                                                                          \n"
	"	/* 0 = xx, 1 = xy, 2 = xz */                                                                                          \n"
	"	/* 3 = yy, 4 = yz, 5 = zz */                                                                                          \n"
	"	float Tij[8];                                                                                                         \n"
	"                                                                                                                         \n"
	"	/* get block and thead indices */                                                                                     \n"
	"	bid = get_group_id( 0);                                                                                               \n"
	"	tid = get_local_id( 0);                                                                                               \n"
	"	gsize = get_num_groups( 0);                                                                                           \n"
	"	bsize = get_local_size( 0);                                                                                           \n"
	"                                                                                                                         \n"
	"	/* this is the induced dipole that we are summing for */                                                              \n"
	"	i = bid*bsize + tid;                                                                                                  \n"
	"                                                                                                                         \n"
	"	/* set the arrays */                                                                                                  \n"
	"	pos = arg_ptr; arg_ptr += N;                                                                                          \n"
	"	estat = arg_ptr; arg_ptr += N;                                                                                        \n"
	"	mu_in = arg_ptr; arg_ptr += N;                                                                                        \n"
	"	mu_out = arg_ptr; arg_ptr += N;                                                                                       \n"
	"	eind_out = arg_ptr;                                                                                                   \n"
	"                                                                                                                         \n"
	"	/* clear the induced field */                                                                                         \n"
	"	eind.x = eind.y = eind.z = eind.w = 0.0f;                                                                             \n"
	"                                                                                                                         \n"
	"	/* locate estat fetch close to pos fetching */                                                                        \n"
	"	estati = estat[i];                                                                                                    \n"
	"	posi = pos[i];                                                                                                        \n"
	"	damp2 = damp*damp;                                                                                                    \n"
	"	damp3 = damp2*damp;                                                                                                   \n"
	"                                                                                                                         \n"
	"	/* for each thread block */                                                                                           \n"
	"	for(j = 0; j < gsize; j++) {                                                                                          \n"
	"                                                                                                                         \n"
	"		/* fill the share mem dipole array for this block */                                                              \n"
	"		smu[tid] = mu_in[j*bsize+tid];                                                                                    \n"
	"		sposj[tid] = pos[j*bsize+tid];                                                                                    \n"
	"                                                                                                                         \n"
	"		/* make sure our shared mem update is complete */                                                                 \n"
	"		barrier( CLK_LOCAL_MEM_FENCE);                                                                                    \n"
	"		/* do work with the shared mem array */                                                                           \n"
	"		#pragma unroll                                                                                                    \n"
	"		for(k = 0; k < BLOCKSIZE; k++) {                                                                                  \n"
	"           barrier( CLK_LOCAL_MEM_FENCE);                                                                                \n"
	"			/* local registers */                                                                                         \n"
	"			mu = smu[k];                                                                                                  \n"
	"			posj = sposj[k];                                                                                              \n"
	"                                                                                                                         \n"
	"			/* START MINIMUM IMAGE */                                                                                     \n"
	"			/* get the particle displacement */                                                                           \n"
	"			dr.x = posi.x - posj.x;                                                                                       \n"
	"			dr.y = posi.y - posj.y;                                                                                       \n"
	"			dr.z = posi.z - posj.z;                                                                                       \n"
	"                                                                                                                         \n"
	"			/* this switch will enforce both i != j and the (MAXATOM-N) null pairs */                                     \n"
	"			sw = (float)(!((int)dr.x) && !((int)dr.y) && !((int)dr.z));                                                   \n"
	"                                                                                                                         \n"
	"			/* matrix multiply with the inverse basis and round */                                                        \n"
	"			img.x = recip_basis_r_a.x*dr.x + recip_basis_r_a.y*dr.y + recip_basis_r_a.z*dr.z;                             \n"
	"			img.y = recip_basis_r_b.x*dr.x + recip_basis_r_b.y*dr.y + recip_basis_r_b.z*dr.z;                             \n"
	"			img.z = recip_basis_r_c.x*dr.x + recip_basis_r_c.y*dr.y + recip_basis_r_c.z*dr.z;                             \n"
	"			img.x = rint(img.x);                                                                                          \n"
	"			img.y = rint(img.y);                                                                                          \n"
	"			img.z = rint(img.z);                                                                                          \n"
	"                                                                                                                         \n"
	"			/* matrix multiply to project back into our basis */                                                          \n"
	"			dri.x = basis_r_a.x*img.x + basis_r_a.y*img.y + basis_r_a.z*img.z;                                            \n"
	"			dri.y = basis_r_b.x*img.x + basis_r_b.y*img.y + basis_r_b.z*img.z;                                            \n"
	"			dri.z = basis_r_c.x*img.x + basis_r_c.y*img.y + basis_r_c.z*img.z;                                            \n"
	"                                                                                                                         \n"
	"			/* now correct the displacement */                                                                            \n"
	"			dri.x = dr.x - dri.x;                                                                                         \n"
	"			dri.y = dr.y - dri.y;                                                                                         \n"
	"			dri.z = dr.z - dri.z;                                                                                         \n"
	"           r2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;                                                                       \n"
	"			ri2 = dri.x*dri.x + dri.y*dri.y + dri.z*dri.z;                                                                \n"
	"           if( !r2 && !ri2) continue;                                                                                    \n"
	"           r = native_sqrt(r2);                                                                                          \n"
	"           ri = native_sqrt(ri2);                                                                                        \n"
	"                                                                                                                         \n"
	"			/* various powers of r that we need */                                                                        \n"
//			"			r2 += sw * MAXFVALUE;                                                                                         \n"
//			"			r  = native_sqrt(r2);                                                                                         \n"
	"           if( isnan(ri) == 0){ r = ri; r2 = ri2;}                                                                       \n"
	"           if( isnan(dri.x) != 0) dri.x = dr.x;                                                                          \n"
	"           if( isnan(dri.y) != 0) dri.y = dr.y;                                                                          \n"
	"           if( isnan(dri.z) != 0) dri.z = dr.z;                                                                          \n"
	"			r3 = r*r*r;                                                                                                   \n"
	"			r5 = r3*r2;                                                                                                   \n"
	"			r3 = native_recip(r3);                                                                                        \n"
	"			r5 = native_recip(r5);                                                                                        \n"
	"			/* END MINIMUM IMAGE */                                                                                       \n"
	"                                                                                                                         \n"
	"			/* damping terms */                                                                                           \n"
	"			expr =  native_exp(-damp*r);                                                                                  \n"
	"			damping_term1 = 1.0f - expr*(0.5f*damp2*r2 + damp*r + 1.0f);                                                  \n"
	"			damping_term2 = 1.0f - expr*(damp3*r*native_divide(r2,6.0f) + 0.5f*damp2*r2 + damp*r + 1.0f);                 \n"
	"                                                                                                                         \n"
	"			/* construct the Tij tensor field, unrolled by hand to avoid conditional on the diagonal terms */             \n"
	"			damping_term1 *= r3;                                                                                          \n"
	"			damping_term2 *= -3.0f*r5;                                                                                    \n"
	"			/* exploit symmetry */                                                                                        \n"
	"			/* 0 = xx, 1 = xy, 2 = xz */                                                                                  \n"
	"			/* 3 = yy, 4 = yz, 5 = zz */                                                                                  \n"
	"			Tij[0] = dri.x*dri.x*damping_term2 + damping_term1;                                                           \n"
	"			Tij[1] = dri.x*dri.y*damping_term2;                                                                           \n"
	"			Tij[2] = dri.x*dri.z*damping_term2;                                                                           \n"
	"			Tij[3] = dri.y*dri.y*damping_term2 + damping_term1;                                                           \n"
	"			Tij[4] = dri.y*dri.z*damping_term2;                                                                           \n"
	"			Tij[5] = dri.z*dri.z*damping_term2 + damping_term1;                                                           \n"
	"                                                                                                                         \n"
	"			/* contract dipole with the tensor */                                                                         \n"
	"			eind.x -= Tij[0]*mu.x + Tij[1]*mu.y + Tij[2]*mu.z;                                                            \n"
	"			eind.y -= Tij[1]*mu.x + Tij[3]*mu.y + Tij[4]*mu.z;                                                            \n"
	"			eind.z -= Tij[2]*mu.x + Tij[4]*mu.y + Tij[5]*mu.z;                                                            \n"
	"                                                                                                                         \n"
	"		} /* end k */                                                                                                     \n"
	"                                                                                                                         \n"
	"	} /* end j */                                                                                                         \n"
	"                                                                                                                         \n"
	"	/* update the ith induced field vector and dipole in global mem */                                                    \n"
	"	mu_out_r.x = estati.w*(estati.x + eind.x);                                                                            \n"
	"	mu_out_r.y = estati.w*(estati.y + eind.y);                                                                            \n"
	"	mu_out_r.z = estati.w*(estati.z + eind.z);                                                                            \n"
	"	mu_out[i] = mu_out_r;                                                                                                 \n"
	"	eind_out[i] = eind;                                                                                                   \n"
	"}		                                                                                                                  \n"
	"__kernel void thole_estat_ocl                                                                                            \n"
	"(                                                                                                                        \n"
	"  __global float4 *arg_ptr,                                                                                              \n"
	"  __global float  *basis,                                                                                                \n"
	"  __global float  *recip_basis,         																				  \n"
	"  __const  uint    N,                                                                                                    \n"
	"  __const  float   pbc_cutoff,                                                                                           \n"
	"  __const  float   gamma,                                                                                                \n"
	"  __const  __global int     *frozen,                                                                                     \n"
	"  __const  __global int     *mol_id                                                                                      \n"
	")                                                                                                                        \n"
	"{                                                                                                                        \n"
	"   __local float4 sposj[ BLOCKSIZE];                                                                                     \n"
	"   __local int sfrozenj[ BLOCKSIZE];                                                                                     \n"
	"   __local int smolidj[ BLOCKSIZE];                                                                                      \n"
	"	int bid, tid, bsize, gsize, mol_id_j, frozen_j;                                                                       \n"
	"        int i, j, k;                                                                                                     \n"
	"	__global float4 *pos, *estat, *mu_in, *mu_out;                                                                        \n"
	"        float4 dr, dri, img;                                                                                             \n"
	"	float4 posi, posj, estati, mu_out_j, mu_out_i;                                                                        \n"
	"	float r, r2, r3, ri, ri2, charge_j, gamma_pol;                                                                        \n"
	"   float damp = DAMP, damp3, damping;                               	/* exponential damping func */                    \n"
	"	float expr, dmp_chrg_r3;                                                                                              \n"
	"	float sw;							/* switch to avoid conditional */                                                 \n"
	"   float4 basis_r_a = (float4)( basis[ 0], basis[ 1], basis[ 2], 0.0f);                                                  \n"
	"   float4 basis_r_b = (float4)( basis[ 3], basis[ 4], basis[ 5], 0.0f);                                                  \n"
	"   float4 basis_r_c = (float4)( basis[ 6], basis[ 7], basis[ 8], 0.0f);                                                  \n"
	"   float4 recip_basis_r_a = (float4)( recip_basis[ 0], recip_basis[ 1], recip_basis[ 2], 0.0f);                          \n"
	"   float4 recip_basis_r_b = (float4)( recip_basis[ 3], recip_basis[ 4], recip_basis[ 5], 0.0f);                          \n"
	"   float4 recip_basis_r_c = (float4)( recip_basis[ 6], recip_basis[ 7], recip_basis[ 8], 0.0f);                          \n"
	"                                                                                                                         \n"
	"	/* get block and thead indices */                                                                                     \n"
	"	bid = get_group_id( 0);                                                                                               \n"
	"	tid = get_local_id( 0);                                                                                               \n"
	"	gsize = get_num_groups( 0);                                                                                           \n"
	"	bsize = get_local_size( 0);                                                                                           \n"
	"                                                                                                                         \n"
	"	/* this is the induced dipole that we are summing for */                                                              \n"
	"	i = bid*bsize + tid;                                                                                                  \n"
	"                                                                                                                         \n"
	"	/* set the arrays */                                                                                                  \n"
	"	pos = arg_ptr; arg_ptr += N;                                                                                          \n"
	"	estat = arg_ptr; arg_ptr += N;                                                                                        \n"
	"	mu_in = arg_ptr; arg_ptr += N;                                                                                        \n"
	"   mu_out = arg_ptr;                                                                                                     \n"
	"                                                                                                                         \n"
	"	/* locate estat fetch close to pos fetching */                                                                        \n"
	"	estati = estat[i];                                                                                                    \n"
	"	posi = pos[i];                                                                                                        \n"
	"   mu_out_i = mu_out[i];                                                                                                 \n"
	"   int mol_id_i = mol_id[i];                                                                                             \n"
	"   int frozen_i = frozen[i];                                                                                             \n"
	"	damp3 = damp*damp*damp;                                                                                               \n"
	"	/* for each thread block */                                                                                           \n"
	"	for(j = 0; j < gsize; j++) {                                                                                          \n"
	"                                                                                                                         \n"
	"		/* fill the share mem dipole array for this block */                                                              \n"
	"		sposj[tid] = pos[j*bsize+tid];                                                                                    \n"
	"       sfrozenj[tid] = frozen[j*bsize+tid];                                                                              \n"
	"       smolidj[tid] = mol_id[j*bsize+tid];                                                                               \n"
	"		/* make sure our shared mem update is complete */                                                                 \n"
	"		barrier( CLK_LOCAL_MEM_FENCE);                                                                                    \n"
	"		/* do work with the shared mem array */                                                                           \n"
	"		#pragma unroll                                                                                                    \n"
	"		for(k = 0; k < BLOCKSIZE; k++) {                                                                                  \n"
	"           barrier( CLK_LOCAL_MEM_FENCE);                                                                                \n"
	"			/* local registers */                                                                                         \n"
	"			posj = sposj[k];                                                                                              \n"
	"           mol_id_j = smolidj[k];                                                                                        \n"
	"           frozen_j = sfrozenj[k];                                                                                       \n"
	"           if( frozen_i && frozen_j) continue;                                                                           \n"
	"           if( mol_id_i == mol_id_j) continue;                                                                           \n"
	"           if( posj.w == 0.0f || posi.w == 0.0f) continue;                                                               \n"
	"           charge_j = posj.w;                                                                                            \n"
	"			/* START MINIMUM IMAGE */                                                                                     \n"
	"			/* get the particle displacement */                                                                           \n"
	"			dr.x = posi.x - posj.x;                                                                                       \n"
	"			dr.y = posi.y - posj.y;                                                                                       \n"
	"			dr.z = posi.z - posj.z;                                                                                       \n"
	"                                                                                                                         \n"
	"			/* this switch will enforce both i != j and the (MAXATOM-N) null pairs */                                     \n"
//			"			sw = (float)(!((int)dr.x) && !((int)dr.y) && !((int)dr.z));                                                   \n"
	"                                                                                                                         \n"
	"			/* matrix multiply with the inverse basis and round */                                                        \n"
	"			img.x = recip_basis_r_a.x*dr.x + recip_basis_r_a.y*dr.y + recip_basis_r_a.z*dr.z;                             \n"
	"			img.y = recip_basis_r_b.x*dr.x + recip_basis_r_b.y*dr.y + recip_basis_r_b.z*dr.z;                             \n"
	"			img.z = recip_basis_r_c.x*dr.x + recip_basis_r_c.y*dr.y + recip_basis_r_c.z*dr.z;                             \n"
	"			img.x = rint(img.x);                                                                                          \n"
	"			img.y = rint(img.y);                                                                                          \n"
	"			img.z = rint(img.z);                                                                                          \n"
	"                                                                                                                         \n"
	"			/* matrix multiply to project back into our basis */                                                          \n"
	"			dri.x = basis_r_a.x*img.x + basis_r_a.y*img.y + basis_r_a.z*img.z;                                            \n"
	"			dri.y = basis_r_b.x*img.x + basis_r_b.y*img.y + basis_r_b.z*img.z;                                            \n"
	"			dri.z = basis_r_c.x*img.x + basis_r_c.y*img.y + basis_r_c.z*img.z;                                            \n"
	"                                                                                                                         \n"
	"			/* now correct the displacement */                                                                            \n"
	"			dri.x = dr.x - dri.x;                                                                                         \n"
	"			dri.y = dr.y - dri.y;                                                                                         \n"
	"			dri.z = dr.z - dri.z;                                                                                         \n"
	"           r2 = dr.x*dr.x + dr.y*dr.y + dr.z*dr.z;                                                                       \n"
	"			ri2 = dri.x*dri.x + dri.y*dri.y + dri.z*dri.z;                                                                \n"
	"           if( !r2 && !ri2) continue;                                                                                    \n"
	"           r = native_sqrt(r2);                                                                                          \n"
	"           ri = native_sqrt(ri2);                                                                                        \n"
	"                                                                                                                         \n"
	"			/* various powers of r that we need */                                                                        \n"
//			"			r2 += sw * MAXFVALUE;                                                                                         \n"
//			"			r  = native_sqrt(r2);                                                                                         \n"
	"           if( isnan(ri) == 0) r = ri;                                                                                   \n"
	"           if(r >= pbc_cutoff || r == 0.0f) continue;                                                                    \n"
	"           if( isnan(dri.x) != 0) dri.x = dr.x;                                                                          \n"
	"           if( isnan(dri.y) != 0) dri.y = dr.y;                                                                          \n"
	"           if( isnan(dri.z) != 0) dri.z = dr.z;                                                                          \n"
	"			r3 = r*r*r;                                                                                                   \n"
	"			/* damping terms */                                                                                           \n"
	"			expr =  native_exp( -native_divide(r3,damp3));                                                                \n"
	"			damping = 1.0f - expr;                                                                                        \n"
	"			/* END MINIMUM IMAGE */                                                                                       \n"
	"                                                                                                                         \n"
	"			r3 = native_recip(r3);                                                                                        \n"
	"           dmp_chrg_r3 = damping * charge_j * r3;                                                                        \n"
	"           estati.x += dmp_chrg_r3 * dri.x;                                                                              \n"
	"           estati.y += dmp_chrg_r3 * dri.y;                                                                              \n"
	"           estati.z += dmp_chrg_r3 * dri.z;                                                                              \n"
	"                                                                                                                         \n"
	"		} /* end k */                                                                                                     \n"
	"                                                                                                                         \n"
	"	} /* end j */                                                                                                         \n"
	"  estat[i] = estati;                                                                                                     \n"
	"  gamma_pol = gamma * estati.w;                                                                                          \n"
	"  mu_in[ i].x = gamma_pol * estati.x;                                                                                    \n"
	"  mu_in[ i].y = gamma_pol * estati.y;                                                                                    \n"
	"  mu_in[ i].z = gamma_pol * estati.z;                                                                                    \n"
	"}		                                                                                                                  \n"
	"__kernel void palmo_echg                                                                                                 \n"
	"(                                                                                                                        \n"
	"  __global float4 *arg_ptr,                                                                                              \n"
	"  __global float4 *echg,                                                                                                 \n"
	"  __global float  *basis,                                                                                                \n"
	"  __global float  *recip_basis,         																				  \n"
	"  __const  uint    N                                                                                                     \n"
	")                                                                                                                        \n"
	"{  																													  \n"
	"   __local float4 sposj[ BLOCKSIZE];                                                                                     \n"
	"   __local float4 smu[ BLOCKSIZE];                                                                                       \n"
	"	int bid, tid, bsize, gsize;                                                                                           \n"
	"        int i, j, k;                                                                                                     \n"
	"	__global float4 *pos, *estat, *mu_in, *mu_out, *eind_out;                                                             \n"
	"        float4 dr, dri, img;                                                                                             \n"
	"	float4 posi, posj, estati, mu, eind;                                                                                  \n"
	"	float r, r2, r3, r5;                                                                                                  \n"
	"        float damp = DAMP, damp2, damp3, damping_term1, damping_term2;	/* exponential damping func */                    \n"
	"	float expr;                                                                                                           \n"
	"	float sw;							/* switch to avoid conditional */                                                 \n"
	"   float4 basis_r_a = (float4)( basis[ 0], basis[ 1], basis[ 2], 0.0f);                                                  \n"
	"   float4 basis_r_b = (float4)( basis[ 3], basis[ 4], basis[ 5], 0.0f);                                                  \n"
	"   float4 basis_r_c = (float4)( basis[ 6], basis[ 7], basis[ 8], 0.0f);                                                  \n"
	"   float4 recip_basis_r_a = (float4)( recip_basis[ 0], recip_basis[ 1], recip_basis[ 2], 0.0f);                          \n"
	"   float4 recip_basis_r_b = (float4)( recip_basis[ 3], recip_basis[ 4], recip_basis[ 5], 0.0f);                          \n"
	"   float4 recip_basis_r_c = (float4)( recip_basis[ 6], recip_basis[ 7], recip_basis[ 8], 0.0f);                          \n"
	"   /* dipole field tensor    */                                                                                          \n"
	"	/* 0 = xx, 1 = xy, 2 = xz */                                                                                          \n"
	"	/* 3 = yy, 4 = yz, 5 = zz */                                                                                          \n"
	"	float Tij[8];                                                                                                         \n"
	"                                                                                                                         \n"
	"	/* get block and thead indices */                                                                                     \n"
	"	bid = get_group_id( 0);                                                                                               \n"
	"	tid = get_local_id( 0);                                                                                               \n"
	"	gsize = get_num_groups( 0);                                                                                           \n"
	"	bsize = get_local_size( 0);                                                                                           \n"
	"                                                                                                                         \n"
	"	/* this is the induced dipole that we are summing for */                                                              \n"
	"	i = get_global_id( 0);                                                                                                \n"
	"	pos = arg_ptr; arg_ptr += N;                                                                                          \n"
	"	estat = arg_ptr; arg_ptr += N;                                                                                        \n"
	"	mu_in = arg_ptr; arg_ptr += N;                                                                                        \n"
	"	mu_out = arg_ptr; arg_ptr += N;                                                                                       \n"
	"	eind_out = arg_ptr;                                                                                                   \n"
	"                                                                                                                         \n"
	"	/* clear the induced field */                                                                                         \n"
	"	eind.x = eind.y = eind.z = eind.w = 0.0f;                                                                             \n"
	"                                                                                                                         \n"
	"	/* locate estat fetch close to pos fetching */                                                                        \n"
	"	estati = estat[i];                                                                                                    \n"
	"	posi = pos[i];                                                                                                        \n"
	"                                                                                                                         \n"
	"	/* for each thread block */                                                                                           \n"
	"	for(j = 0; j < gsize; j++) {                                                                                          \n"
	"                                                                                                                         \n"
	"		/* fill the share mem dipole array for this block */                                                              \n"
	"       int idx = mad24(j,bsize,tid);                                                                                     \n"
	"		smu[tid] = mu_in[idx];                                                                                            \n"
	"		sposj[tid] = pos[idx];                                                                                            \n"
	"                                                                                                                         \n"
	"		/* make sure our shared mem update is complete */                                                                 \n"
	"		barrier( CLK_LOCAL_MEM_FENCE);                                                                                    \n"
	"		/* do work with the shared mem array */    																		  \n"
	"       #pragma unroll                                                          										  \n"
	"		for(k = 0; k < BLOCKSIZE; k++) {                                                                                  \n"
	"           barrier( CLK_LOCAL_MEM_FENCE);                                                                                \n"
	"			/* local registers */                                                                                         \n"
	"			mu = smu[k];                                                                                                  \n"
	"			posj = sposj[k];                                                                                              \n"
	"                                                                                                                         \n"
	"			/* START MINIMUM IMAGE */                                                                                     \n"
	"			/* get the particle displacement */                                                                           \n"
	"			dr.x = posi.x - posj.x;                                                                                       \n"
	"			dr.y = posi.y - posj.y;                                                                                       \n"
	"			dr.z = posi.z - posj.z;                                                                                       \n"
	"			/* this switch will enforce both i != j and the (MAXATOM-N) null pairs */                                     \n"
	"			sw = (float)(!((int)dr.x) && !((int)dr.y) && !((int)dr.z));                                                   \n"
	"                                                                                                                         \n"
	"			/* matrix multiply with the inverse basis and round */                                                        \n"
	"			img.x = recip_basis_r_a.x*dr.x + recip_basis_r_a.y*dr.y + recip_basis_r_a.z*dr.z;                             \n"
	"			img.y = recip_basis_r_b.x*dr.x + recip_basis_r_b.y*dr.y + recip_basis_r_b.z*dr.z;                             \n"
	"			img.z = recip_basis_r_c.x*dr.x + recip_basis_r_c.y*dr.y + recip_basis_r_c.z*dr.z;                             \n"
	"			img.x = rint(img.x);                                                                                          \n"
	"			img.y = rint(img.y);                                                                                          \n"
	"			img.z = rint(img.z);                                                                                          \n"
	"                                                                                                                         \n"
	"			/* matrix multiply to project back into our basis */                                                          \n"
	"			dri.x = basis_r_a.x*img.x + basis_r_a.y*img.y + basis_r_a.z*img.z;                                            \n"
	"			dri.y = basis_r_b.x*img.x + basis_r_b.y*img.y + basis_r_b.z*img.z;                                            \n"
	"			dri.z = basis_r_c.x*img.x + basis_r_c.y*img.y + basis_r_c.z*img.z;                                            \n"
	"                                                                                                                         \n"
	"			/* now correct the displacement */                                                                            \n"
	"			dri.x = dr.x - dri.x;                                                                                         \n"
	"			dri.y = dr.y - dri.y;                                                                                         \n"
	"			dri.z = dr.z - dri.z;                                                                                         \n"
	"			r2 = dri.x*dri.x + dri.y*dri.y + dri.z*dri.z;                                                                 \n"
	"                                                                                                                         \n"
	"			/* various powers of r that we need */                                                                        \n"
	"			r2 += sw * MAXFVALUE;                                                                                         \n"
	"			r  = native_sqrt(r2);                                                                                         \n"
	"			r3 = r2*r;                                                                                                    \n"
	"			r5 = r3*r2;                                                                                                   \n"
	"			r3 = native_recip(r3);                                                                                        \n"
	"			r5 = native_recip(r5);                                                                                        \n"
	"			/* END MINIMUM IMAGE */                                                                                       \n"
	"                                                                                                                         \n"
	"			/* damping terms */                                                                                           \n"
	"			damp2 = damp*damp;                                                                                            \n"
	"			damp3 = damp2*damp;                                                                                           \n"
	"			expr =  native_exp(-damp*r);                                                                                  \n"
	"			damping_term1 = 1.0f - expr*(0.5f*damp2*r2 + damp*r + 1.0f);                                                  \n"
	"			damping_term2 = 1.0f - expr*(damp3*r*native_divide(r2,6.0f) + 0.5f*damp2*r2 + damp*r + 1.0f);                 \n"
	"                                                                                                                         \n"
	"			/* construct the Tij tensor field, unrolled by hand to avoid conditional on the diagonal terms */             \n"
	"			damping_term1 *= r3;                                                                                          \n"
	"			damping_term2 *= -3.0f*r5;                                                                                    \n"
	"			/* exploit symmetry */                                                                                        \n"
	"			/* 0 = xx, 1 = xy, 2 = xz */                                                                                  \n"
	"			/* 3 = yy, 4 = yz, 5 = zz */                                                                                  \n"
	"			Tij[0] = dri.x*dri.x*damping_term2 + damping_term1;                                                           \n"
	"			Tij[1] = dri.x*dri.y*damping_term2;                                                                           \n"
	"			Tij[2] = dri.x*dri.z*damping_term2;                                                                           \n"
	"			Tij[3] = dri.y*dri.y*damping_term2 + damping_term1;                                                           \n"
	"			Tij[4] = dri.y*dri.z*damping_term2;                                                                           \n"
	"			Tij[5] = dri.z*dri.z*damping_term2 + damping_term1;                                                           \n"
	"                                                                                                                         \n"
	"			/* contract dipole with the tensor */                                                                         \n"
	"			eind.x -= Tij[0]*mu.x + Tij[1]*mu.y + Tij[2]*mu.z;                                                            \n"
	"			eind.y -= Tij[1]*mu.x + Tij[3]*mu.y + Tij[4]*mu.z;                                                            \n"
	"			eind.z -= Tij[2]*mu.x + Tij[4]*mu.y + Tij[5]*mu.z;                                                            \n"
	"                                                                                                                         \n"
	"		} /* end k */                                                                                                     \n"
	"                                                                                                                         \n"
	"	} /* end j */                                                                                                         \n"
	"                                                                                                                         \n"
	"	echg[i] = eind;                                                                                                       \n"
	"}   	                                                                                                                  \n"
	"__kernel void PotentialReduction                                                                                         \n"
	"(                                                                                                                        \n"
	"  __global float4 *arg_ptr,                                                                                              \n"
	"  __const  uint    N,                                                                                                    \n"
	"  __global float  *OUTPUT                                                                                                \n"
	")                                                                                                                        \n"
	"{                                                                                                                        \n"
	"  __local float SHARED[ BLOCKSIZE];                                                                                      \n"
	"  // perform first level of reduction,                                                                                   \n"
	"  // reading from global memory, writing to shared memory                                                                \n"
	"  const uint thread_id = get_local_id( 0);                                                                               \n"
	"  const uint block_size = get_local_size( 0);                                                                            \n"
	"  const uint grid_size = block_size * 2 * get_num_groups( 0);                                                            \n"
	"  uint index = get_group_id( 0) * get_local_size( 0) * 2 + get_local_id( 0);                                             \n"
	"  SHARED[ thread_id] = 0;                                                                                                \n"
	"  __global float4 *estat, *mu_out;                                                                                       \n"
	"  arg_ptr += N;                                                                                                          \n"
	"  estat = arg_ptr; arg_ptr += 2*N;                                                                                       \n"
	"  mu_out = arg_ptr;                                                                                                      \n"
	"                                                                                                                         \n"
	"  // we reduce multiple elements per thread.  The number is determined by the                                            \n"
	"  // number of active thread blocks (via gridDim).  More blocks will result                                              \n"
	"  // in a larger grid_size and therefore fewer elements per thread                                                       \n"
	"  while( index < N)                                                                                                      \n"
	"  {                                                                                                                      \n"
	"	SHARED[ thread_id] += estat[ index].x * mu_out[ index].x;                                                             \n"
	"	SHARED[ thread_id] += estat[ index].y * mu_out[ index].y;                                                             \n"
	"	SHARED[ thread_id] += estat[ index].z * mu_out[ index].z;                                                             \n"
	"	// ensure we don't read out of bounds                                                                                 \n"
	"	if( index + block_size < N)                                                                                           \n"
	"	{                                                                                                                     \n"
	"	  SHARED[ thread_id] += estat[ index + block_size].x * mu_out[ index + block_size].x;                                 \n"
	"	  SHARED[ thread_id] += estat[ index + block_size].y * mu_out[ index + block_size].y;                                 \n"
	"	  SHARED[ thread_id] += estat[ index + block_size].z * mu_out[ index + block_size].z;                                 \n"
	"	}                                                                                                                     \n"
	"	index += grid_size;                                                                                                   \n"
	"  }                                                                                                                      \n"
	"                                                                                                                         \n"
	"  barrier( CLK_LOCAL_MEM_FENCE);                                                                                         \n"
	"                                                                                                                         \n"
	"  // do reduction in shared mem                                                                                          \n"
	"  for( uint offset = block_size / 2; offset > 0; offset = offset / 2)                                                    \n"
	"  {                                                                                                                      \n"
	"	if( thread_id < offset)                                                                                               \n"
	"	{                                                                                                                     \n"
	"	  SHARED[ thread_id] += SHARED[ thread_id + offset];                                                                  \n"
	"	}                                                                                                                     \n"
	"	barrier( CLK_LOCAL_MEM_FENCE);                                                                                        \n"
	"  }                                                                                                                      \n"
	"                                                                                                                         \n"
	"  // write result for this block to global mem                                                                           \n"
	"  if( thread_id == 0)                                                                                                    \n"
	"  {                                                                                                                      \n"
	"	OUTPUT[ get_group_id( 0)] = SHARED[ 0];                                                                               \n"
	"  }                                                                                                                      \n"
	"}                                                                                                                        \n"
	;

//	fprintf( stderr, "after kernel_code\n");
	unsigned int nr_platforms;
	unsigned int nr_devices;
	cl_int error = CL_SUCCESS;
	// platform
	error = clGetPlatformIDs( 0, NULL, &nr_platforms);
	CheckErr( error, "clGetPlatformIDs");
	ocl->platforms = calloc( nr_platforms, sizeof( cl_platform_id));
	error = clGetPlatformIDs( nr_platforms, ocl->platforms, &nr_platforms);
	CheckErr( error, "clGetPlatformIDs");

//	fprintf( stderr, "found %i platforms\n", nr_platforms);
	int ctr;
	for( ctr = 0; ctr < nr_platforms; ctr++)
	{
//		fprintf( stderr, "checking platform %i\n\n", ctr);
		error = clGetDeviceIDs( ocl->platforms[ ctr], CL_DEVICE_TYPE_GPU, 0, NULL, &nr_devices);
//		CheckErr( error, "clGetDeviceIDs");
//		fprintf( stderr, "number devices found of type GPU: %i\n\n", nr_devices);
		if( error != CL_SUCCESS && error != CL_DEVICE_NOT_FOUND)
		{
//			fprintf( stderr, "in if\n");
			CheckErr( error, "clGetDeviceIDs");
		}
		else if( error == CL_SUCCESS && nr_devices > 0)
		{
//			fprintf( stderr, "in else if\n");
			ocl->device_id = calloc( nr_devices, sizeof( cl_device_id));
			error = clGetDeviceIDs( ocl->platforms[ ctr], CL_DEVICE_TYPE_GPU, nr_devices, &ocl->device_id[ 0], NULL);
			CheckErr( error, "clGetDeviceIDs");
			break;
		}
	}

	// Check to see if we found at least one GPU device, otherwise return
	if( ocl->device_id == NULL)
	{
		fprintf( stderr, "No GPU device found\n");
		exit(-1);
	}

	// context
	cl_context_properties properties[] = {CL_CONTEXT_PLATFORM, ( cl_context_properties)( ocl->platforms[ctr]), 0};
	ocl->context = clCreateContextFromType( properties, CL_DEVICE_TYPE_GPU, NULL, NULL, &error);
	CheckErr( error, "clCreateContextFromType");

	size_t data;
	size_t kernel_length;

	ocl->queue = clCreateCommandQueue( ocl->context, ocl->device_id[ 0], 0, &error);
	CheckErr( error, "clCreateCommandQueue");
	size_t lengths[ 1];
	lengths[ 0] = strlen( const_kernel_code);

	ocl->program = clCreateProgramWithSource( ocl->context, 1, &const_kernel_code, lengths, &error);
	CheckErr( error, "clCreateProgramWithSource");

	error = clBuildProgram( ocl->program, 1, &ocl->device_id[ 0], NULL, NULL, NULL);
	CheckErr( error, "clBuildProgram");
//	cl_build_status build_status;
//	error = clGetProgramBuildInfo( ocl->program, ocl->device_id[ 0], CL_PROGRAM_BUILD_STATUS, sizeof( cl_build_status), &build_status, NULL);
//    char *build_log;
//    size_t ret_val_size;
//    error = clGetProgramBuildInfo(ocl->program, ocl->device_id[ 0], CL_PROGRAM_BUILD_LOG, 0, NULL, &ret_val_size);
//    CheckErr(error, "clGetProgramBuildInfo");
//    build_log = malloc(ret_val_size+1);
//    error = clGetProgramBuildInfo(ocl->program, ocl->device_id[ 0], CL_PROGRAM_BUILD_LOG, ret_val_size, build_log, NULL);
//    CheckErr(error, "clGetProgramBuildInfo");
//
//    // to be carefully, terminate with \0
//    // there's no information in the reference whether the string is 0 terminated or not
//    build_log[ret_val_size] = '\0';
//
//    fprintf(stderr, "%s\n\n", build_log );

	ocl->kernel = clCreateKernel( ocl->program, "thole_iterative_ocl", &error);
	CheckErr( error, "clCreateKernel");
	ocl->palmo_echg = clCreateKernel( ocl->program, "palmo_echg", &error);
	CheckErr( error, "clCreateKernel");
	ocl->thole_estat = clCreateKernel( ocl->program, "thole_estat_ocl", &error);
	CheckErr( error, "clCreateKernel");
	ocl->potential_reduction = clCreateKernel( ocl->program, "PotentialReduction", &error);
	CheckErr( error, "clCreateKernel");

	return ocl;
}
#endif

